% A Basic Mortensen-Pissarides Type of Search Model with SGU First-Order
% Local Approximation on Hagedorn-Manovskii Model
% Written by: Orhan Torul
% Istanbul, Turkey, May 2017
%%
%%
%Housekeeping
clc;
clear all;
close all;

%%
%%
%Parameters

global beta delta xi kappa psi nu rhoz sigmaz sigma;
global zss;

beta=0.99;
delta=0.10; 
sigma=1;
xi=0.40;
kappa=0.50;
psi=0.25;
nu=0.40;
rhoz=0.95;
sigmaz=0.007;

%%
%%
%Steady-State
%c  : consumption
%u  : unemployment rate (n = 1-u)
%v  : vacancy postings
%q  : rate at which job vacations are filled
%p  : rate at which jobs are filled
%w  : wage
%R  : rate of return
%chi: unemployment benefit/home production determined as a ratio of w (0.4 or 0.95)

zss=0;

%y=[c; u ;v ;q ;p ;R ;w ;chi]
y0=[0.96 ;   0.37  ;  0.06  ;  1.02  ;  0.18  ;  1.01 ;   0.97  ;  0.92];
%defining optimization options
options=optimset('Display','iter','MaxIter',10000, 'MaxFunEvals', 10000,'tolfun',10e-25,'tolx',10e-25,'FunValCheck','on','PlotFcns',@optimplotx);
%calling for optimization
[out,fval,exitflag] = fsolve(@steadyhm,y0,options); 
if exitflag>0;
    disp('Steady-State Derivation Successful!');
end

css=out(1);
uss=out(2);
vss=out(3);
thetass=out(3)/out(2);
qss=out(4);
pss=out(5);
Rss=out(6);
wss=out(7);
chiss=out(8);


%Derivation of the Preliminaries of SGU Matrices

syms c u v p q chi w R z Xi;
syms c_1 u_1 v_1 q_1 p_1 z_1 w_1 Xi_1 R_1 chi_1; 

f1=-psi/(q*(1-delta))+Xi_1*((exp(z_1))-w_1+psi/q_1); 
%f1 is the job posting condition
f2=u*p-v*q;
%f2 is the equivalence of the match function definitions
f3=1-R*beta*((c_1/c)^(-sigma));
%f3 is the interest rate equation
f4=-w+nu*((exp(z))+psi*(v/u))+(1-nu)*chi;
%f4 is the wage as a result of wage bargaining
f5=-chi+0.95*w;
%f5 is the Shimer condition
f6=-(1-u_1)+(1-delta)*(1-u)+v*q;
%f6 is the evolution of the labor force 
f7=kappa*(u^xi)*(v^(1-xi))-v*q;
%f7 is the match function relation with q&v (and implicitly u&p)
f8=exp(z)*(1-u)+u*chi-c-psi*v;
%f8 is the aggregate resource constraint
f9=Xi_1-beta*((c_1/c)^(-sigma));
%f9 is the definition of discount rate Xi
f10=z_1-(1-rhoz)*zss-rhoz*z-randn(1,1)*sigmaz;
%f10 is the evolution of Solow residual/technology shock

F=[f1; f2; f3; f4; f5; f6; f7; f8; f9 ; f10];
%F is the vectorized form of all equilibrium conditions

y = [c; v; p; q; chi; w; R; Xi]; 
%y is the co-state vector as of time t
y_1= [c_1; v_1; p_1; q_1; chi_1; w_1; R_1; Xi_1];
%y_1 is the co-state vector as of time t+1
x = [u; z]; 
%x is the state vector as of time t
x_1= [u_1; z_1];
%x is the state vector as of time t+1


f_y=jacobian(F,y); 
f_y_1=jacobian(F,y_1);
f_x=jacobian(F,x);
f_x_1=jacobian(F,x_1);
%Above are the Jacobians to be used for SGU alghorithm.

%Setting the variables to their long-run (steady-state) values
c=css;
u=uss;
v=vss;
q=qss;
p=pss;
R=Rss;
w=wss;
chi=chiss;
z=zss;
Xiss=beta;

Xi=beta;
c_1=css; 
u_1=uss; 
v_1=vss;
q_1=qss;
p_1=pss;
z_1=zss;
w_1=wss;
chi_1=chiss;
R_1=Rss;
Xi_1=beta;

f_y=eval(f_y);
f_y_1=eval(f_y_1);
f_x=eval(f_x);
f_x_1=eval(f_x_1);
%Above are the steady-state values of the Jacobians

global ga ha;
syms a11g a12g a21g a22g a31g a32g a41g a42g a51g a52g a61g a62g a71g a72g a81g a82g;
syms a11h a12h a21h a22h;

g=[a11g a12g; a21g a22g; a31g a32g; a41g a42g; a51g a52g; a61g a62g; a71g a72g; a81g a82g];
h=[a11h a12h; a21h a22h];

G =f_y_1*g*h + f_y*g + f_x_1*h + f_x;

Gvector=G(:);

Ginitial=[0.2;0.2;0.3;0.3;0.1;0.1;0.5;0.5;0.4;0.4;0.2;0.2;0.3;0.3;0.1;0.1;0.5;0.5;0.4;0.4;];
options=optimset('Display','iter','MaxIter',10000, 'MaxFunEvals', 10000,'tolfun',10e-25,'tolx',10e-25,'FunValCheck','on','PlotFcns',@optimplotx);

[matrices,fval,exitflag] = fsolve(@(varhm) sguhm(varhm,f_y_1,f_y,f_x_1,f_x),Ginitial,options); 
if exitflag>0;
    disp('SGU Derivation Successful!');
end
%%
%%
%Simulations
Z=zeros(100,100);
load shocks;
TFP=SHOCKS;
randn('seed', 2703);
TFP(:,70)=randn(1,100);

%Reminder: The order for the ha matrix is as follows= [u;z]; 
for i=1:100;
    for t=1:99;
        Z(t+1,i)=(1-rhoz)*zss+rhoz*Z(t,i)+sigmaz*TFP(t+1,i);
    end;
end;

U=zeros(100,100);
for j=1:100;
    U(1,j)=uss;
    for t=1:99;
        U(t+1,j)=uss+(U(t,j)-uss)*ha(1,1)+(Z(t,j)-zss)*ha(1,2);
    end;
end;


C=zeros(100,100);
V=zeros(100,100);
P=zeros(100,100);
Q=zeros(100,100);
CHI=zeros(100,100);
W=zeros(100,100);
RR=zeros(100,100);
XI=zeros(100,100);

%Reminder: The order for the ga matrix is as follows= [c; v; p; q; chi; w; R; Xi]; 
for j=1:100;
    for t=1:100;
        C(t,j)=css+(U(t,j)-uss)*ga(1,1)+(Z(t,j)-zss)*ga(1,2);
        V(t,j)=vss+(U(t,j)-uss)*ga(2,1)+(Z(t,j)-zss)*ga(2,2);
        P(t,j)=pss+(U(t,j)-uss)*ga(3,1)+(Z(t,j)-zss)*ga(3,2);
        Q(t,j)=qss+(U(t,j)-uss)*ga(4,1)+(Z(t,j)-zss)*ga(4,2);
        CHI(t,j)=chiss+(U(t,j)-uss)*ga(5,1)+(Z(t,j)-zss)*ga(5,2);
        W(t,j)=wss+(U(t,j)-uss)*ga(6,1)+(Z(t,j)-zss)*ga(6,2);
        RR(t,j)=Rss+(U(t,j)-uss)*ga(7,1)+(Z(t,j)-zss)*ga(7,2);
        XI(t,j)=Xiss+(U(t,j)-uss)*ga(8,1)+(Z(t,j)-zss)*ga(8,2);
    end
end
THETA=V./U;
N=1-U;
OUTPUT=exp(Z).*N;


%HP Filter
[Uf]=hpfilter(log(U),1600);
Ur=U-Uf;
[Zf]=hpfilter(Z,1600);
Zr=Z-Zf;
[Cf]=hpfilter(log(C),1600);
Cr=C-Cf;
[Vf]=hpfilter(log(V),1600);
Vr=V-Vf;
[Pf]=hpfilter(log(P),1600);
Pr=P-Pf;
[Qf]=hpfilter(log(Q),1600);
Qr=Q-Qf;
[CHIf]=hpfilter(log(CHI),1600);
CHIr=CHI-CHIf;
[Wf]=hpfilter(log(W),1600);
Wr=W-Wf;
[RRf]=hpfilter(log(RR),1600);
RRr=RR-RRf;
[THETAf]=hpfilter(log(THETA),1600);
THETAr=THETA-THETAf;
[Nf]=hpfilter(log(N),1600);
Nr=N-Nf;
[OUTPUTf]=hpfilter(log(OUTPUT),1600);
OUTPUTr=OUTPUT-OUTPUTf;

% [Uf]=hpfilter(log(U),10000);
% Ur=U-Uf;
% [Zf]=hpfilter(Z,10000);
% Zr=Z-Zf;
% [Cf]=hpfilter(log(C),10000);
% Cr=C-Cf;
% [Vf]=hpfilter(log(V),10000);
% Vr=V-Vf;
% [Pf]=hpfilter(log(P),10000);
% Pr=P-Pf;
% [Qf]=hpfilter(log(Q),10000);
% Qr=Q-Qf;
% [CHIf]=hpfilter(log(CHI),10000);
% CHIr=CHI-CHIf;
% [Wf]=hpfilter(log(W),10000);
% Wr=W-Wf;
% [RRf]=hpfilter(log(RR),10000);
% RRr=RR-RRf;
% [THETAf]=hpfilter(log(THETA),10000);
% THETAr=THETA-THETAf;
% [Nf]=hpfilter(log(N),10000);
% Nr=N-Nf;
% [OUTPUTf]=hpfilter(log(OUTPUT),10000);
% OUTPUTr=OUTPUT-OUTPUTf;



%Descriptive Statistics
meanU=zeros(100,1);
meanZ=zeros(100,1);
meanC=zeros(100,1);
meanV=zeros(100,1);
meanP=zeros(100,1);
meanQ=zeros(100,1);
meanCHI=zeros(100,1);
meanW=zeros(100,1);
meanRR=zeros(100,1);
meanTHETA=zeros(100,1);
meanN=zeros(100,1);
meanOUTPUT=zeros(100,1);

stdevU=zeros(100,1);
stdevZ=zeros(100,1);
stdevC=zeros(100,1);
stdevV=zeros(100,1);
stdevP=zeros(100,1);
stdevQ=zeros(100,1);
stdevCHI=zeros(100,1);
stdevW=zeros(100,1);
stdevRR=zeros(100,1);
stdevTHETA=zeros(100,1);
stdevN=zeros(100,1);
stdevOUTPUT=zeros(100,1);

corrZ_U=zeros(100,1);
corrZ_Z=zeros(100,1);
corrZ_C=zeros(100,1);
corrZ_V=zeros(100,1);
corrZ_P=zeros(100,1);
corrZ_Q=zeros(100,1);
corrZ_CHI=zeros(100,1);
corrZ_W=zeros(100,1);
corrZ_RR=zeros(100,1);
corrZ_THETA=zeros(100,1);
corrZ_N=zeros(100,1);
corrZ_OUTPUT=zeros(100,1);

corrOUTPUT_U=zeros(100,1);
corrOUTPUT_Z=zeros(100,1);
corrOUTPUT_C=zeros(100,1);
corrOUTPUT_V=zeros(100,1);
corrOUTPUT_P=zeros(100,1);
corrOUTPUT_Q=zeros(100,1);
corrOUTPUT_CHI=zeros(100,1);
corrOUTPUT_W=zeros(100,1);
corrOUTPUT_RR=zeros(100,1);
corrOUTPUT_THETA=zeros(100,1);
corrOUTPUT_N=zeros(100,1);
corrOUTPUT_OUTPUT=zeros(100,1);



corrU_U=zeros(100,1);
corrU_Z=zeros(100,1);
corrU_C=zeros(100,1);
corrU_V=zeros(100,1);
corrU_P=zeros(100,1);
corrU_Q=zeros(100,1);
corrU_CHI=zeros(100,1);
corrU_W=zeros(100,1);
corrU_RR=zeros(100,1);
corrU_THETA=zeros(100,1);
corrU_N=zeros(100,1);
corrU_OUTPUT=zeros(100,1);

corrV_U=zeros(100,1);
corrV_Z=zeros(100,1);
corrV_C=zeros(100,1);
corrV_V=zeros(100,1);
corrV_P=zeros(100,1);
corrV_Q=zeros(100,1);
corrV_CHI=zeros(100,1);
corrV_W=zeros(100,1);
corrV_RR=zeros(100,1);
corrV_THETA=zeros(100,1);
corrV_N=zeros(100,1);
corrV_OUTPUT=zeros(100,1);

corrP_U=zeros(100,1);
corrP_Z=zeros(100,1);
corrP_C=zeros(100,1);
corrP_V=zeros(100,1);
corrP_P=zeros(100,1);
corrP_Q=zeros(100,1);
corrP_CHI=zeros(100,1);
corrP_W=zeros(100,1);
corrP_RR=zeros(100,1);
corrP_THETA=zeros(100,1);
corrP_N=zeros(100,1);
corrP_OUTPUT=zeros(100,1);

corrQ_U=zeros(100,1);
corrQ_Z=zeros(100,1);
corrQ_C=zeros(100,1);
corrQ_V=zeros(100,1);
corrQ_P=zeros(100,1);
corrQ_Q=zeros(100,1);
corrQ_CHI=zeros(100,1);
corrQ_W=zeros(100,1);
corrQ_RR=zeros(100,1);
corrQ_THETA=zeros(100,1);
corrQ_N=zeros(100,1);
corrQ_OUTPUT=zeros(100,1);

corrTHETA_U=zeros(100,1);
corrTHETA_Z=zeros(100,1);
corrTHETA_C=zeros(100,1);
corrTHETA_V=zeros(100,1);
corrTHETA_P=zeros(100,1);
corrTHETA_Q=zeros(100,1);
corrTHETA_CHI=zeros(100,1);
corrTHETA_W=zeros(100,1);
corrTHETA_RR=zeros(100,1);
corrTHETA_THETA=zeros(100,1);
corrTHETA_N=zeros(100,1);
corrTHETA_OUTPUT=zeros(100,1);

coefofvarU=zeros(100,1);
coefofvarZ=zeros(100,1);
coefofvarC=zeros(100,1);
coefofvarV=zeros(100,1);
coefofvarP=zeros(100,1);
coefofvarQ=zeros(100,1);
coefofvarCHI=zeros(100,1);
coefofvarW=zeros(100,1);
coefofvarRR=zeros(100,1);
coefofvarTHETA=zeros(100,1);
coefofvarN=zeros(100,1);
coefofvarOUTPUT=zeros(100,1);

autocorrU=zeros(100,2);
autocorrZ=zeros(100,2);
autocorrC=zeros(100,2);
autocorrV=zeros(100,2);
autocorrP=zeros(100,2);
autocorrQ=zeros(100,2);
autocorrCHI=zeros(100,2);
autocorrW=zeros(100,2);
autocorrRR=zeros(100,2);
autocorrTHETA=zeros(100,2);
autocorrN=zeros(100,2);
autocorrOUTPUT=zeros(100,2);


for i=1:100;
    meanU(i,1)=mean(Ur(:,i));
    meanZ(i,1)=mean(Zr(:,i));
    meanC(i,1)=mean(Cr(:,i));
    meanV(i,1)=mean(Vr(:,i));
    meanP(i,1)=mean(Pr(:,i));
    meanQ(i,1)=mean(Qr(:,i));
    meanCHI(i,1)=mean(CHIr(:,i));
    meanW(i,1)=mean(Wr(:,i));
    meanRR(i,1)=mean(RRr(:,i));
    meanTHETA(i,1)=mean(THETAr(:,i));
    meanN(i,1)=mean(Nr(:,i));
    meanOUTPUT(i,1)=mean(OUTPUTr(:,i));
end


for i=1:100;
    stdevU(i,1)=std(Ur(:,i));
    stdevZ(i,1)=std(Zr(:,i));
    stdevC(i,1)=std(Cr(:,i));
    stdevV(i,1)=std(Vr(:,i));
    stdevP(i,1)=std(Pr(:,i));
    stdevQ(i,1)=std(Qr(:,i));
    stdevCHI(i,1)=std(CHIr(:,i));
    stdevW(i,1)=std(Wr(:,i));
    stdevRR(i,1)=std(RRr(:,i));
    stdevTHETA(i,1)=std(THETAr(:,i));
    stdevN(i,1)=std(Nr(:,i));
    stdevOUTPUT(i,1)=std(OUTPUTr(:,i));
end


for i=1:100;
    coefofvarU(i,1)=std(Ur(:,i))/mean(Ur(:,i));
    coefofvarZ(i,1)=std(Zr(:,i))/mean(Zr(:,i));
    coefofvarC(i,1)=std(Cr(:,i))/mean(Cr(:,i));
    coefofvarV(i,1)=std(Vr(:,i))/mean(Vr(:,i));
    coefofvarP(i,1)=std(Pr(:,i))/mean(Pr(:,i));
    coefofvarQ(i,1)=std(Qr(:,i))/mean(Qr(:,i));
    coefofvarCHI(i,1)=std(CHIr(:,i))/mean(CHIr(:,i));
    coefofvarW(i,1)=std(Wr(:,i))/mean(Wr(:,i));
    coefofvarRR(i,1)=std(RRr(:,i))/mean(RRr(:,i));
    coefofvarTHETA(i,1)=std(THETAr(:,i))/mean(THETAr(:,i));
    coefofvarN(i,1)=std(Nr(:,i))/mean(Nr(:,i));
    coefofvarOUTPUT(i,1)=std(OUTPUTr(:,i))/mean(OUTPUTr(:,i));
end

for i=1:100;
    corrOUTPUT_U(i,1)=corr(Ur(:,i),OUTPUTr(:,i));
    corrOUTPUT_Z(i,1)=corr(Zr(:,i),OUTPUTr(:,i));
    corrOUTPUT_C(i,1)=corr(Cr(:,i),OUTPUTr(:,i));
    corrOUTPUT_V(i,1)=corr(Vr(:,i),OUTPUTr(:,i));
    corrOUTPUT_P(i,1)=corr(Pr(:,i),OUTPUTr(:,i));
    corrOUTPUT_Q(i,1)=corr(Qr(:,i),OUTPUTr(:,i));
    corrOUTPUT_CHI(i,1)=corr(CHIr(:,i),OUTPUTr(:,i));
    corrOUTPUT_W(i,1)=corr(Wr(:,i),OUTPUTr(:,i));
    corrOUTPUT_RR(i,1)=corr(RRr(:,i),OUTPUTr(:,i));
    corrOUTPUT_THETA(i,1)=corr(THETAr(:,i),OUTPUTr(:,i));
    corrOUTPUT_N(i,1)=corr(Nr(:,i),OUTPUTr(:,i));
    corrOUTPUT_OUTPUT(i,1)=corr(OUTPUTr(:,i),OUTPUTr(:,i));
end


for i=1:100;
    corrZ_U(i,1)=corr(Ur(:,i),Zr(:,i));
    corrZ_Z(i,1)=corr(Zr(:,i),Zr(:,i));
    corrZ_C(i,1)=corr(Cr(:,i),Zr(:,i));
    corrZ_V(i,1)=corr(Vr(:,i),Zr(:,i));
    corrZ_P(i,1)=corr(Pr(:,i),Zr(:,i));
    corrZ_Q(i,1)=corr(Qr(:,i),Zr(:,i));
    corrZ_CHI(i,1)=corr(CHIr(:,i),Zr(:,i));
    corrZ_W(i,1)=corr(Wr(:,i),Zr(:,i));
    corrZ_RR(i,1)=corr(RRr(:,i),Zr(:,i));
    corrZ_THETA(i,1)=corr(THETAr(:,i),Zr(:,i));
    corrZ_N(i,1)=corr(Nr(:,i),Zr(:,i));
    corrZ_OUTPUT(i,1)=corr(OUTPUTr(:,i),Zr(:,i));
end

for i=1:100;
    corrU_U(i,1)=corr(Ur(:,i),Ur(:,i));
    corrU_Z(i,1)=corr(Zr(:,i),Ur(:,i));
    corrU_C(i,1)=corr(Cr(:,i),Ur(:,i));
    corrU_V(i,1)=corr(Vr(:,i),Ur(:,i));
    corrU_P(i,1)=corr(Pr(:,i),Ur(:,i));
    corrU_Q(i,1)=corr(Qr(:,i),Ur(:,i));
    corrU_CHI(i,1)=corr(CHIr(:,i),Ur(:,i));
    corrU_W(i,1)=corr(Wr(:,i),Ur(:,i));
    corrU_RR(i,1)=corr(RRr(:,i),Ur(:,i));
    corrU_THETA(i,1)=corr(THETAr(:,i),Ur(:,i));
    corrU_N(i,1)=corr(Nr(:,i),Ur(:,i));
    corrU_OUTPUT(i,1)=corr(OUTPUTr(:,i),Ur(:,i));
end

for i=1:100;
    corrV_U(i,1)=corr(Ur(:,i),Vr(:,i));
    corrV_Z(i,1)=corr(Zr(:,i),Vr(:,i));
    corrV_C(i,1)=corr(Cr(:,i),Vr(:,i));
    corrV_V(i,1)=corr(Vr(:,i),Vr(:,i));
    corrV_P(i,1)=corr(Pr(:,i),Vr(:,i));
    corrV_Q(i,1)=corr(Qr(:,i),Vr(:,i));
    corrV_CHI(i,1)=corr(CHIr(:,i),Vr(:,i));
    corrV_W(i,1)=corr(Wr(:,i),Vr(:,i));
    corrV_RR(i,1)=corr(RRr(:,i),Vr(:,i));
    corrV_THETA(i,1)=corr(THETAr(:,i),Vr(:,i));
    corrV_N(i,1)=corr(Nr(:,i),Vr(:,i));
    corrV_OUTPUT(i,1)=corr(OUTPUTr(:,i),Vr(:,i));
end

for i=1:100;
    corrP_U(i,1)=corr(Ur(:,i),Pr(:,i));
    corrP_Z(i,1)=corr(Zr(:,i),Pr(:,i));
    corrP_C(i,1)=corr(Cr(:,i),Pr(:,i));
    corrP_V(i,1)=corr(Vr(:,i),Pr(:,i));
    corrP_P(i,1)=corr(Pr(:,i),Pr(:,i));
    corrP_Q(i,1)=corr(Qr(:,i),Pr(:,i));
    corrP_CHI(i,1)=corr(CHIr(:,i),Pr(:,i));
    corrP_W(i,1)=corr(Wr(:,i),Pr(:,i));
    corrP_RR(i,1)=corr(RRr(:,i),Pr(:,i));
    corrP_THETA(i,1)=corr(THETAr(:,i),Pr(:,i));
    corrP_N(i,1)=corr(Nr(:,i),Pr(:,i));
    corrP_OUTPUT(i,1)=corr(OUTPUTr(:,i),Pr(:,i));
end

for i=1:100;
    corrQ_U(i,1)=corr(Ur(:,i),Qr(:,i));
    corrQ_Z(i,1)=corr(Zr(:,i),Qr(:,i));
    corrQ_C(i,1)=corr(Cr(:,i),Qr(:,i));
    corrQ_V(i,1)=corr(Vr(:,i),Qr(:,i));
    corrQ_P(i,1)=corr(Pr(:,i),Qr(:,i));
    corrQ_Q(i,1)=corr(Qr(:,i),Qr(:,i));
    corrQ_CHI(i,1)=corr(CHIr(:,i),Qr(:,i));
    corrQ_W(i,1)=corr(Wr(:,i),Qr(:,i));
    corrQ_RR(i,1)=corr(RRr(:,i),Qr(:,i));
    corrQ_THETA(i,1)=corr(THETAr(:,i),Qr(:,i));
    corrQ_N(i,1)=corr(Nr(:,i),Qr(:,i));
    corrQ_OUTPUT(i,1)=corr(OUTPUTr(:,i),Qr(:,i));
end

for i=1:100;
    corrTHETA_U(i,1)=corr(Ur(:,i),THETAr(:,i));
    corrTHETA_Z(i,1)=corr(Zr(:,i),THETAr(:,i));
    corrTHETA_C(i,1)=corr(Cr(:,i),THETAr(:,i));
    corrTHETA_V(i,1)=corr(Vr(:,i),THETAr(:,i));
    corrTHETA_P(i,1)=corr(Pr(:,i),THETAr(:,i));
    corrTHETA_Q(i,1)=corr(Qr(:,i),THETAr(:,i));
    corrTHETA_CHI(i,1)=corr(CHIr(:,i),THETAr(:,i));
    corrTHETA_W(i,1)=corr(Wr(:,i),THETAr(:,i));
    corrTHETA_RR(i,1)=corr(RRr(:,i),THETAr(:,i));
    corrTHETA_THETA(i,1)=corr(THETAr(:,i),THETAr(:,i));
    corrTHETA_N(i,1)=corr(Nr(:,i),THETAr(:,i));
    corrTHETA_OUTPUT(i,1)=corr(OUTPUTr(:,i),THETAr(:,i));
end

for i=1:100;
    [autocorrU(i,:)]=autocorr(Ur(:,i),1);
    [autocorrZ(i,:)]=autocorr(Zr(:,i),1);
    [autocorrC(i,:)]=autocorr(Cr(:,i),1);
    [autocorrV(i,:)]=autocorr(Vr(:,i),1);
    [autocorrP(i,:)]=autocorr(Pr(:,i),1);
    [autocorrQ(i,:)]=autocorr(Qr(:,i),1);
    [autocorrCHI(i,:)]=autocorr(CHIr(:,i),1);
    [autocorrW(i,:)]=autocorr(Wr(:,i),1);
    [autocorrRR(i,:)]=autocorr(RRr(:,i),1);
    [autocorrTHETA(i,:)]=autocorr(THETAr(:,i),1);
    [autocorrN(i,:)]=autocorr(Nr(:,i),1);
    [autocorrOUTPUT(i,:)]=autocorr(OUTPUTr(:,i),1);
end

MU=mean(meanU);
MZ=mean(meanZ);
MC=mean(meanC);
MV=mean(meanV);
MP=mean(meanP);
MQ=mean(meanQ);
MCHI=mean(meanCHI);
MW=mean(meanW);
MRR=mean(meanRR);
MTHETA=mean(meanTHETA);
MN=mean(meanN);
MOUTPUT=mean(meanOUTPUT);

STDU=mean(stdevU);
STDZ=mean(stdevZ);
STDC=mean(stdevC);
STDV=mean(stdevV);
STDP=mean(stdevP);
STDQ=mean(stdevQ);
STDCHI=mean(stdevCHI);
STDW=mean(stdevW);
STDRR=mean(stdevRR);
STDTHETA=mean(stdevTHETA);
STDN=mean(stdevN);
STDOUTPUT=mean(stdevOUTPUT);

CORROUTPUTU=mean(corrOUTPUT_U);
CORROUTPUTZ=mean(corrOUTPUT_Z);
CORROUTPUTC=mean(corrOUTPUT_C);
CORROUTPUTV=mean(corrOUTPUT_V);
CORROUTPUTP=mean(corrOUTPUT_P);
CORROUTPUTQ=mean(corrOUTPUT_Q);
CORROUTPUTCHI=mean(corrOUTPUT_CHI);
CORROUTPUTW=mean(corrOUTPUT_W);
CORROUTPUTRR=mean(corrOUTPUT_RR);
CORROUTPUTTHETA=mean(corrOUTPUT_THETA);
CORROUTPUTN=mean(corrOUTPUT_N);
CORROUTPUTOUTPUT=mean(corrOUTPUT_OUTPUT);

CORRZU=mean(corrZ_U);
CORRZZ=mean(corrZ_Z);
CORRZC=mean(corrZ_C);
CORRZV=mean(corrZ_V);
CORRZP=mean(corrZ_P);
CORRZQ=mean(corrZ_Q);
CORRZCHI=mean(corrZ_CHI);
CORRZW=mean(corrZ_W);
CORRZRR=mean(corrZ_RR);
CORRZTHETA=mean(corrZ_THETA);
CORRZN=mean(corrZ_N);
CORRZOUTPUT=mean(corrZ_OUTPUT);


CORRUU=mean(corrU_U);
CORRUZ=mean(corrU_Z);
CORRUC=mean(corrU_C);
CORRUV=mean(corrU_V);
CORRUP=mean(corrU_P);
CORRUQ=mean(corrU_Q);
CORRUCHI=mean(corrU_CHI);
CORRUW=mean(corrU_W);
CORRURR=mean(corrU_RR);
CORRUTHETA=mean(corrU_THETA);
CORRUN=mean(corrU_N);
CORRUOUTPUT=mean(corrU_OUTPUT);

CORRVU=mean(corrV_U);
CORRVZ=mean(corrV_Z);
CORRVC=mean(corrV_C);
CORRVV=mean(corrV_V);
CORRVP=mean(corrV_P);
CORRVQ=mean(corrV_Q);
CORRVCHI=mean(corrV_CHI);
CORRVW=mean(corrV_W);
CORRVRR=mean(corrV_RR);
CORRVTHETA=mean(corrV_THETA);
CORRVN=mean(corrV_N);
CORRVOUTPUT=mean(corrV_OUTPUT);

CORRPU=mean(corrP_U);
CORRPZ=mean(corrP_Z);
CORRPC=mean(corrP_C);
CORRPV=mean(corrP_V);
CORRPP=mean(corrP_P);
CORRPQ=mean(corrP_Q);
CORRPCHI=mean(corrP_CHI);
CORRPW=mean(corrP_W);
CORRPRR=mean(corrP_RR);
CORRPTHETA=mean(corrP_THETA);
CORRPN=mean(corrP_N);
CORRPOUTPUT=mean(corrP_OUTPUT);

CORRQU=mean(corrQ_U);
CORRQZ=mean(corrQ_Z);
CORRQC=mean(corrQ_C);
CORRQV=mean(corrQ_V);
CORRQP=mean(corrQ_P);
CORRQQ=mean(corrQ_Q);
CORRQCHI=mean(corrQ_CHI);
CORRQW=mean(corrQ_W);
CORRQRR=mean(corrQ_RR);
CORRQTHETA=mean(corrQ_THETA);
CORRQN=mean(corrQ_N);
CORRQOUTPUT=mean(corrQ_OUTPUT);

CORRTHETAU=mean(corrTHETA_U);
CORRTHETAZ=mean(corrTHETA_Z);
CORRTHETAC=mean(corrTHETA_C);
CORRTHETAV=mean(corrTHETA_V);
CORRTHETAP=mean(corrTHETA_P);
CORRTHETAQ=mean(corrTHETA_Q);
CORRTHETACHI=mean(corrTHETA_CHI);
CORRTHETAW=mean(corrTHETA_W);
CORRTHETARR=mean(corrTHETA_RR);
CORRTHETATHETA=mean(corrTHETA_THETA);
CORRTHETAN=mean(corrTHETA_N);
CORRTHETAOUTPUT=mean(corrTHETA_OUTPUT);


CVARU=mean(coefofvarU);
CVARZ=mean(coefofvarZ);
CVARC=mean(coefofvarC);
CVARV=mean(coefofvarV);
CVARP=mean(coefofvarP);
CVARQ=mean(coefofvarQ);
CVARCHI=mean(coefofvarCHI);
CVARW=mean(coefofvarW);
CVARRR=mean(coefofvarRR);
CVARTHETA=mean(coefofvarTHETA);
CVARN=mean(coefofvarN);
CVAROUTPUT=mean(coefofvarOUTPUT);

ACORRU=mean(autocorrU(:,2));
ACORRZ=mean(autocorrZ(:,2));
ACORRC=mean(autocorrC(:,2));
ACORRV=mean(autocorrV(:,2));
ACORRP=mean(autocorrP(:,2));
ACORRQ=mean(autocorrQ(:,2));
ACORRCHI=mean(autocorrCHI(:,2));
ACORRW=mean(autocorrW(:,2));
ACORRRR=mean(autocorrRR(:,2));
ACORRTHETA=mean(autocorrTHETA(:,2));
ACORRN=mean(autocorrN(:,2));
ACORROUTPUT=mean(autocorrOUTPUT(:,2));

RELU=STDU/STDOUTPUT;
RELZ=STDZ/STDOUTPUT;
RELC=STDC/STDOUTPUT;
RELV=STDV/STDOUTPUT;
RELP=STDP/STDOUTPUT;
RELQ=STDQ/STDOUTPUT;
RELCHI=STDCHI/STDOUTPUT;
RELW=STDW/STDOUTPUT;
RELRR=STDRR/STDOUTPUT;
RELTHETA=STDTHETA/STDOUTPUT;
RELN=STDN/STDOUTPUT;
RELOUTPUT=STDOUTPUT/STDOUTPUT;

RELUZ=STDU/STDZ;
RELZZ=STDZ/STDZ;
RELCZ=STDC/STDZ;
RELVZ=STDV/STDZ;
RELPZ=STDP/STDZ;
RELQZ=STDQ/STDZ;
RELCHIZ=STDCHI/STDZ;
RELWZ=STDW/STDZ;
RELRRZ=STDRR/STDZ;
RELTHETAZ=STDTHETA/STDZ;
RELNZ=STDN/STDZ;
RELOUTPUTZ=STDOUTPUT/STDZ;


%Importing Appendix to LaTeX
Meanmatrix=[MU MZ MC MV MP MQ MW MRR MTHETA MOUTPUT];
Stdmatrix=[STDU STDZ STDC STDV STDP STDQ STDW STDRR STDTHETA  STDOUTPUT];
Corrtoymatrix=[CORROUTPUTU CORROUTPUTZ CORROUTPUTC CORROUTPUTV CORROUTPUTP CORROUTPUTQ  CORROUTPUTW CORROUTPUTRR CORROUTPUTTHETA  CORROUTPUTOUTPUT];
Corrtoymatrixz=[CORRZU CORRZZ CORRZC CORRZV CORRZP CORRZQ  CORRZW CORRZRR CORRZTHETA  CORRZOUTPUT];
%Coefofvarmatrix=[CVARU CVARZ CVARC CVARV CVARP CVARQ  CVARW CVARRR CVARTHETA CVAROUTPUT ];
Relstdmatrix=[RELU RELZ RELC RELV RELP RELQ RELW RELRR RELTHETA RELOUTPUT];
Relstdmatrixz=[RELUZ RELZZ RELCZ RELVZ RELPZ RELQZ RELWZ RELRRZ RELTHETAZ RELOUTPUTZ];
Autocorrmatrix=[ACORRU ACORRZ ACORRC ACORRV ACORRP ACORRQ  ACORRW ACORRRR ACORRTHETA ACORROUTPUT]; 
    
TableDescriptive=[Stdmatrix ; Relstdmatrixz; Relstdmatrix; Corrtoymatrixz; Corrtoymatrix ; Autocorrmatrix];
rowlabeldesc={'Std. Deviation', 'Rel. Std. ($e^z$)','Rel. Std. ($y$)','Corr.($e^z$)', 'Corr.($y$)', 'Autocorrelation'};
columnlabeldesc = {'$u$', '$e^z$' , '$c$' ,'$v$', '$p$', '$q$',  '$w$' ,'$R$' , '$\theta$', '$y$'};
matrix2latex(TableDescriptive, 'tableDescriptive_hm.tex','rowLabels', rowlabeldesc, 'columnLabels', columnlabeldesc,'format', '%-6.4f', 'size', 'footnotesize');

Corrumatrix=[CORRUU CORRUV CORRUTHETA CORRUP CORRUQ];
Corrvmatrix=[CORRVU CORRVV CORRVTHETA CORRVP CORRVQ];
Corrthetamatrix=[CORRTHETAU CORRTHETAV CORRTHETATHETA CORRTHETAP CORRTHETAQ];
Corrpmatrix=[CORRPU CORRPV CORRPTHETA CORRPP CORRPQ];
Corrqmatrix=[CORRQU CORRQV CORRQTHETA CORRQP CORRQQ];

CORRMATRIX=[Corrumatrix; Corrvmatrix; Corrthetamatrix; Corrpmatrix; Corrqmatrix];
rowlabelcorr={'$u$', '$v$','$\theta$','$p$','$q$'};
columnlabelcorr={'$u$', '$v$','$\theta$','$p$','$q$'};
matrix2latex(CORRMATRIX, 'corrmatrix_hm.tex','rowLabels', rowlabelcorr, 'columnLabels', columnlabelcorr,'format', '%-6.4f', 'size', 'footnotesize');


nss=1-uss;
outputss=exp(zss)*nss;
Valuesmatrix=[uss zss css vss pss qss chiss wss Rss thetass nss outputss];
ssmatrix=Valuesmatrix;
rowlabeldescss={'SS'};
columnlabeldescss = {'$u$', '$z$' , '$c$' ,'$v$', '$p$', '$q$','$\chi$', '$w$' ,'$R$' , '$\theta$', '$n$', '$y$'};
matrix2latex(ssmatrix, 'tabless_hm.tex','rowLabels', rowlabeldescss, 'columnLabels', columnlabeldescss,'format', '%-6.4f', 'size', 'footnotesize');


%%
%Simulations
Z=zeros(100,100);
load shocks;
TFP=SHOCKS;
randn('seed', 2703);
TFP(:,70)=randn(1,100);

%Reminder: The order for the ha matrix is as follows= [u;z]; 
for i=1:100;
    for t=1:99;
        Z(t+1,i)=0;
    end;
end;

for i=1:100;
    
        Z(1,i)=2*sigmaz;;
 
end;


U=zeros(100,100);

for j=1:100;
    U(1,j)=uss;
    for t=1:99;
        U(t+1,j)=uss+(U(t,j)-uss)*ha(1,1)+(Z(t,j)-zss)*ha(1,2);
    end;
end;



C=zeros(100,100);
V=zeros(100,100);
P=zeros(100,100);
Q=zeros(100,100);
CHI=zeros(100,100);
W=zeros(100,100);
RR=zeros(100,100);
XI=zeros(100,100);

%Reminder: The order for the ga matrix is as follows= [c; v; p; q; chi; w; R; Xi]; 
for j=1:100;
    for t=1:100;
        C(t,j)=css+(U(t,j)-uss)*ga(1,1)+(Z(t,j)-zss)*ga(1,2);
        V(t,j)=vss+(U(t,j)-uss)*ga(2,1)+(Z(t,j)-zss)*ga(2,2);
        P(t,j)=pss+(U(t,j)-uss)*ga(3,1)+(Z(t,j)-zss)*ga(3,2);
        Q(t,j)=qss+(U(t,j)-uss)*ga(4,1)+(Z(t,j)-zss)*ga(4,2);
        CHI(t,j)=chiss+(U(t,j)-uss)*ga(5,1)+(Z(t,j)-zss)*ga(5,2);
        W(t,j)=wss+(U(t,j)-uss)*ga(6,1)+(Z(t,j)-zss)*ga(6,2);
        RR(t,j)=Rss+(U(t,j)-uss)*ga(7,1)+(Z(t,j)-zss)*ga(7,2);
        XI(t,j)=Xiss+(U(t,j)-uss)*ga(8,1)+(Z(t,j)-zss)*ga(8,2);
    end
end
THETA=V./U;
N=1-U;
OUTPUT=exp(Z).*N;



fig=figure;
subplot(1,1,1);
plot(U(2:end,1), '--r', 'LineWidth', 2.5);
ylim([0.357 0.367])
%title('Unemployment', 'fontweight', 'bold','FontSize',14);
xlabel('Period','FontSize',14);
ylabel('Unemployment','FontSize',14);
xt = get(gca, 'XTick');
set(gca, 'FontSize', 14);
grid;
print(fig,'u_hm','-dpng')


fig=figure;
subplot(1,1,1);
plot(V(2:end,1), '-.b', 'LineWidth', 2.5);
ylim([0.06095 0.0625])
%title('Vacancies', 'fontweight', 'bold','FontSize',14);
xlabel('Period','FontSize',14);
ylabel('Vacancy','FontSize',14);
xt = get(gca, 'XTick');
set(gca, 'FontSize', 14);
grid;
print(fig,'v_hm','-dpng')

